*	See Gibbs estimation of microstructure models: Teaching notes
	November, 2006
	Joel Hasbrouck
____________________________________________________________________________________________________

	RollGibbs 01.sas

____________________________________________________________________________________________________;

*	Modify the following statement to set current directory to location of program;
x 'cd c:\Active\Mkt Microstructure Book\Programs Fall 2006';
options nodate nocenter nonumber ps=70 ls=120  sasautos=('./Sas Macros','.') mprint; 
libname SASData '.\Datasets';
libname this '.';

/*ods pdf file='./RollGibbs 01.pdf' style=serifPrinter;*/

%let Infinity=1e30;
%let eps=1e-30;

*___________________________________________________________________________________________________

	Generate some price data
____________________________________________________________________________________________________;
%let nObs=30;
%let sdu = .01;
%let c = .001;
%let c = 0.01;
title "Simulation with sdu=&sdu and c=&c";
data sim;
	m = 0;
	do t=1 to &nObs;
		m = m + &sdu*normal(345699);
		q = 1-2*(uniform(246)>.5);
		p = m + &c*q;
		keep q p;
		output;
	end;
	run;

%let nSweeps=10000;

*___________________________________________________________________________________________________

	Simulate parameters and q's
____________________________________________________________________________________________________;
/*%RollGibbs(sim, nSweeps=&nSweeps, cStart=&c, varuStart=%sysevalf(&sdu*&sdu));*/
*	To estimate the model with no non-negativity restriction on c, use the following;
%RollGibbs(sim, nSweeps=&nSweeps, cStart=&c, varuStart=%sysevalf(&sdu*&sdu),cLower=-&Infinity);


*___________________________________________________________________________________________________

	Analyze the simulated estimates
____________________________________________________________________________________________________;
%let pDrop=0.2;							*	Proportion of sweeps to drop (burn in);
%let nDrop = %sysevalf(&pDrop*&nSweeps);	*	Number of sweeps to drop;
data parmOut2;
	set parmOut;
	if sweep>&nDrop then output;
	run;
proc means data=parmOut2 fw=6;
	run;
proc transpose data=parmOut2 out=parmOut2T (rename=(col1=value));
	by sweep;
	run;
proc sort data=parmOut2T;
	by _name_ sweep;
	run;
proc kde data=parmOut2T;
	univar value / levels out=myOut;
	by _name_;
	run; 
goptions ftext=complex htext=2 colors=(black) hby=0;
symbol1 v=none i=sm w=2;
proc gplot data=myOut;
	plot density*value=1;
	title "Simulation with sdu=&sdu and c=&c: Posterior for #byval1";
	by _name_;
	run;
	quit;
proc kde data=parmOut2;
	bivar  sdu c / levels out=myOut;
	run;
proc gcontour data=myOut;
	plot value1*value2=density;
	label value1='sdu' value2='c';
	title "Joint posterior with sdu=&sdu and c=&c";
	run;
	quit;

/*data qOut2;*/
/*	set qOut;*/
/*	if sweep>&nDrop then output;*/
/*	run;*/
/*proc transpose data=qOut2 out=qt prefix=q;*/
/*	var q;*/
/*	by sweep;*/
/*	run;*/
/*proc means data=qt (drop=sweep) noprint;*/
/*	output out=qtMeans (drop=_type_ _freq_) mean= ;*/
/*	run;*/
/*proc transpose data=qtMeans out=qMeans (rename=(col1=q));*/
/*	run;*/
/*data qMeans;*/
/*	set qMeans;*/
/*	t = _n_;*/
/*	keep t q;*/
/*	run;*/
/*proc print data=qMeans;*/
/*	run;*/
